home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 2.iso / STUTTGART / LANG / SCHEME / GNU / SCM4E1 / !Scm / slib / randinex < prev    next >
Text File  |  1994-03-06  |  4KB  |  99 lines
  1. ;;;; Pseudo-Random inexact real numbers for scheme.
  2. ;;; Copyright (C) 1991, 1993 Aubrey Jaffer.
  3. ;;; New sphere and normal functions from: Harald Hanche-Olsen
  4.  
  5. ;Permission to copy this software, to redistribute it, and to use it
  6. ;for any purpose is granted, subject to the following restrictions and
  7. ;understandings.
  8.  
  9. ;1.  Any copy made of this software must include this copyright notice
  10. ;in full.
  11.  
  12. ;2.  I have made no warrantee or representation that the operation of
  13. ;this software will be error-free, and I am under no obligation to
  14. ;provide any services, by way of maintenance, update, or otherwise.
  15.  
  16. ;3.  In conjunction with products arising from the use of this
  17. ;material, there shall be no use of my name in any advertising,
  18. ;promotional, or sales literature without prior written consent in
  19. ;each case.
  20.  
  21. ;This file is loaded by random.scm if inexact numbers are supported by
  22. ;the implementation.
  23.  
  24. (define random:float-radix
  25.   (+ 1 (exact->inexact random:MASK)))
  26.  
  27. ;;; This determines how many chunks will be neccessary to completely
  28. ;;; fill up an inexact real.
  29. (define (random:size-float l x)
  30.   (if (= 1.0 (+ 1 x))
  31.       l
  32.       (random:size-float (+ l 1) (/ x random:float-radix))))
  33. (define random:chunks/float (random:size-float 0 1.0))
  34.  
  35. (define (random:uniform-chunk n state)
  36.   (if (= 1 n)
  37.       (/ (exact->inexact (random:chunk state))
  38.      random:float-radix)
  39.       (/ (+ (random:uniform-chunk (- n 1) state)
  40.         (exact->inexact (random:chunk state)))
  41.      random:float-radix)))
  42.  
  43. ;;; Generate an inexact real between 0 and 1.
  44. (define (random:uniform state)
  45.   (random:uniform-chunk random:chunks/float state))
  46.  
  47. ;;; If x and y are independent standard normal variables, then with
  48. ;;; x=r*cos(t), y=r*sin(t), we find that t is uniformly distributed
  49. ;;; over [0,2*pi] and the cumulative distribution of r is
  50. ;;; 1-exp(-r^2/2).  This latter means that u=exp(-r^2/2) is uniformly
  51. ;;; distributed on [0,1], so r=sqrt(-2 log u) can be used to generate r.
  52.  
  53. (define (random:normal-vector! vect . args)
  54.   (let ((state (if (null? args) *random-state* (car args)))
  55.     (sum2 0))
  56.     (let ((do! (lambda (k x)
  57.          (vector-set! vect k x)
  58.          (set! sum2 (+ sum2 (* x x))))))
  59.       (do ((n (- (vector-length vect) 1) (- n 2)))
  60.       ((negative? n) sum2)
  61.     (let ((t (* 6.28318530717958 (random:uniform state)))
  62.           (r (sqrt (* -2 (log (random:uniform state))))))
  63.       (do! n (* r (cos t)))
  64.       (if (positive? n) (do! (- n 1) (* r (sin t)))))))))
  65.  
  66. (define random:normal
  67.   (let ((vect (make-vector 1)))
  68.     (lambda args 
  69.       (apply random:normal-vector! vect args)
  70.       (vector-ref vect 0))))
  71.  
  72. ;;; For the uniform distibution on the hollow sphere, pick a normal
  73. ;;; family and scale.
  74.  
  75. (define (random:hollow-sphere! vect . args)
  76.   (let ((ms (sqrt (apply random:normal-vector! vect args))))
  77.     (do ((n (- (vector-length vect) 1) (- n 1)))
  78.     ((negative? n))
  79.       (vector-set! vect n (/ (vector-ref vect n) ms)))))
  80.  
  81. ;;; For the uniform distribution on the solid sphere, note that in
  82. ;;; this distribution the length r of the vector has cumulative
  83. ;;; distribution r^n; i.e., u=r^n is uniform [0,1], so r kan be
  84. ;;; generated as r=u^(1/n).
  85.  
  86. (define (random:solid-sphere! vect . args)
  87.   (apply random:hollow-sphere! vect args)
  88.   (let ((r (expt (random:uniform (if (null? args) *random-state* (car args)))
  89.          (/ (vector-length vect)))))
  90.     (do ((n (- (vector-length vect) 1) (- n 1)))
  91.     ((negative? n))
  92.       (vector-set! vect n (* r (vector-ref vect n))))))
  93.  
  94. (define (random:exp . args)
  95.   (let ((state (if (null? args) *random-state* (car args))))
  96.     (- (log (random:uniform state)))))
  97.  
  98. (require 'random)
  99.